Tarea 2

Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

• Nicolás Acevedo
• Diego Navarrete

Cuerpo Docente:

• Profesor: Felipe Bravo M.
• Auxiliares: Martín Paredes, María José Zambrano
• Ayudantes: Nicolas Cabello, Kevin Iturra, Mariana Vásquez

Índice:

1. Objetivo
2. Instrucciones
3. Referencias
4. Elaboración de Código

Objetivo

Bienvenid@s a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

• La tarea se realiza en grupos de máximo 2 personas. Pero no existe problema si usted desea hacerla de forma individual.
• La entrega es a través de u-cursos a más tardar el día estipulado en la misma plataforma. A las tareas atrasadas se les descontará un punto por día.
• El formato de entrega es este mismo Rmarkdown y un html con la tarea desarrollada. Por favor compruebe que todas las celdas han sido ejecutadas en el archivo html.
• Al momento de la revisión tu código será ejecutado. Por favor verifica que tu entrega no tenga errores de compilación.
• No serán revisadas tareas desarrolladas en Python.
• Está PROHIBIDO la copia o compartir las respuestas entre integrantes de diferentes grupos.
• Pueden realizar consultas de la tarea a través de U-cursos y/o del canal de Discord del curso.

Referencias:

Slides de las clases:

• Introduction to Statistical Inference
• Design of Experiments & Hypothesis Testing

Enlaces a videos de las clases:

• Introduction to Statistical Inference: video1 video2 video3 video4 video5 video6 video7 video8
• Design of Experiments & Hypothesis Testing: video1 video2 video3 video4 video5 video6

Documentación:

• Boostrap en r
• ggplot2

Elaboración de código

Usted deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

Pregunta 1: Estimadores.

Esta pregunta tiene como objetivo comprobar a través de gráficos las características que poseen los estimadores.Por favor responda de forma separada las siguientes preguntas:

• En clases se vio que el estimador \(\hat{p}_{n} = \frac{1}{n} \displaystyle{\sum_{i=1}^{n}}X_{i}\) es un estimador consistente para \(X_{i}\) Bernoulli de tasa \(p\), verifique esto númericamente para una distribución Bernoulli de \(p=0.5\), es decir grafique como se ve \(\hat{p}_{n}\) para valores de \(n\) y compárelo con el valor verdadero.

• Sabemos que no todos los estimadores insesgados son consistentes. Considere el estimador \(T_{n} = \hat{p}_{n} + \epsilon_{n}\) donde \(\epsilon_{n} \sim \mathcal{N}(0,1)\) es posible verificar que \(T_{n}\) es insesgado pero no es consistente, para ver esto repita lo que realizo en el punto anterior y estudie lo que sucede. ¿Porque cree que no es consistente el estimador?, Justifique su respuesta.

Respuesta: En primera instancia, realizamos 50.000 experimentos aleatorios para una distribución de Bernoulli. Luego, graficamos como se va viendo el estimador \(\hat{p}_{n}\) mientras se va avanzando en la cantidad de experimentos.

n_max <- 50000
p <- 0.5

# Vectores del estimador y de la distribución de Bernoulli
estimator <- vector(length = n_max)
bernoulli <- rbinom(n_max, size = 1, prob = p)

# Para cada valor del vector, se promedia con los anteriores
# y se agregan al estimador
for (k in 1:n_max) {
  estimator[k] <- sum(bernoulli[1:k]) / k
}


# Gráfico de los experimentos
plot_est <- plot_ly(
  x = 1:n_max,
  y = estimator,
  mode = "lines",
  type = "scatter",
  name = "Experimentos",
  line = list(color = "blue")
) %>%
add_trace(
  x = 1:n_max,
  y = 0.5,
  name = "p = 0.5",
  line = list(color = "red")
) %>%
layout(
  title = "Estimador Pn para X ~ Bernoulli(p) con p = 0.5",
  xaxis = list(title = sprintf("Pasos (%d)", n_max)),
  yaxis = list(title = "Probabilidad")
)
plot_est

Podemos notar que a medida que vamos avanzando en la cantidad de experimentos, el estimador se acerca cada vez más al valor teórico de \(p = 0.5\), lo cual verifica que el estimador es consistente.

Por otra parte, realizamos el mismo procedimiento para el estimador \(T_{n}\).

n_max <- 50000

estimator <- vector(length = n_max)
bernoulli <- rbinom(n_max, size = 1, prob = 0.5)
noise <- rnorm(n_max)
T <- vector(length = n_max)

for (k in 1:n_max) {
  estimator[k] <- sum(bernoulli[1:k]) / k
  T[k] = estimator[k] + noise[k]
}

plot_est <- plot_ly(
  x = 1:n_max,
  y = T,
  mode = "lines",
  type = "scatter",
  name = "Experimentos",
  line = list(color = "blue")
) %>%
add_trace(
  x = 1:n_max,
  y = 0.5,
  name = "p = 0.5",
  line = list(color = "red")
) %>%
layout(
  title = "Estimador Tn para X ~ Bernoulli(p) con p = 0.5",
  xaxis = list(title = sprintf("Pasos (%d)", n_max)),
  yaxis = list(title = "Probabilidad")
)

plot_est

Se puede apreciar que el estimador \(T_n\) no es consistente ya que no converge al valor esperado de \(p = 0.5\). La inconsistencia de este estimador se debe a que al estimador \(\hat{p}_n\) se le agrega una variable normal en cada iteración, un “ruido”.

plot_est <- plot_ly(
  x = 1:n_max,
  y = noise,
  mode = "lines",
  type = "scatter",
  name = "Ruido",
  line = list(color = "blue")
) %>%
layout(
  title = "Ruido e ~ N(0, 1)",
  xaxis = list(title = sprintf("Pasos (%d)", n_max)),
  yaxis = list(title = "Probabilidad")
)

plot_est

Este ruido es el que no permite que el estimador \(T_n\) converga al parámetro \(p\), ya que hace que el estimador “oscile” aleatoriamente, evitando que converga y por lo tanto, haciendo de \(T_n\) un estimador insesgado pero inconsistente.

---

Pregunta 2: Intervalo de Confianza

El objetivo de esta pregunta es visualizar los intervalos de confianza en datos simulados de una población, para visualizar la incertidumbre que presenta una estimación. Para esto, ustedes deberán generar datos de una distribución exponencial, la cual deberán considerar como los datos de la población. En base a los datos generados, estimen la distribución de la media de la población a través de la sampling distribution de la media. Notar que el valor obtenido en cada muestra les entregará un estimador de la media, o sea, para cada valor podremos calcular un intervalo de confianza. Hecho esto, calculen el intervalo de confianza del \(95\%\) para cada una de las medias estimadas, utilizando la función de cuartil vista en clases.

Para la elaboración de esta parte de la tarea, se recomienda realizar el experimento con la siguiente secuencialidad:

• Obtener la media de la población (en este caso la exponencial).
• Realizar una sampling distribution con un tamaño de muestra igual a \(30\) sobre los datos generados de la población. Repita la obtención de la media un número elevado de veces (recomendación \(5000\) veces).
• Calcular el intervalo de confianza del \(95\%\) para cada uno de las medias obtenidas en las iteraciones.
• De acuerdo a los valores obtenidos (medias e intervalos de confianza), grafique cada una de las medias obtenidas en conjunto a sus intervalos de confianza. Aquí debe notar que, si el intervalo de confianza contiene a la media de la población, este se considerará como parte del intervalo de confianza del \(95\%\), haga un conteo de cuantos valores están contenidos en este intervalo.
• Señalar en el plot de intervalo de confianza los valores que están dentro y fuera del intervalo. Comente que visualiza de los intervalos de confianza obtenidos, ¿existe incertidumbre?.
• En base al conteo realizado en el punto anterior, observe cómo se comporta la proporción de intervalos de confianza, ¿es acaso este igual al \(95\%\) teórico usado para calcularlo?.

Notar: Responder cada una de las preguntas señaladas en esta sección.

Hints:

• Para realizar la sampling distribution podría serle útil el comando sample().
• Puede ser útil generar un dataframe para graficar con ggplot2.

Gráfico esperado para intervalos de confianza

Del gráfico es posible observar que la línea punteada es la media de la población y los puntos de colores son las estimaciones con sus respectivos intervalos de confianza. Notar que para el plot no se utilizaron las 5000 veces, se recomienda utilizar 100 valores para visualizar bien el fenómeno.

---

Respuesta:

#########100 Simulaciones##############
# Paso 1: Generar datos de una distribución exponencial (población)
tamano_poblacion <- 10000
poblacion <- rexp(tamano_poblacion, rate = 0.5)
media_poblacion <- mean(poblacion)
print(media_poblacion)

## [1] 1.979917

# Paso 2: Realizar una sampling distribution de la media
tamano_muestra <- 30
num_simulaciones <- 100
medias_muestrales <- c()

for (i in 1:num_simulaciones) {
  muestra_indices <- sample(poblacion, size = tamano_muestra, replace = FALSE)
  medias_muestrales[i] <- mean(muestra_indices)
}


# Paso 3: Calcular el intervalo de confianza del 95%
alpha <- 0.05
sigma <- 2
se <- sigma/sqrt(tamano_muestra)
error <- qnorm(1 - alpha/2) * se
left <- medias_muestrales - error
right <- medias_muestrales + error


# Paso 4: Crear un dataframe para ggplot2
data <- data.frame(
  Value = medias_muestrales,
  Label = c(1:num_simulaciones),
  LowerCI = left,
  UpperCI = right,
  central = media_poblacion  # Media de los datos (mu)
)
valor_central = data$central

# Paso 5: Graficar las medias obtenidas junto con sus intervalos de confianza
# Crea el gráfico de intervalo de confianza en vertical con segmentos de error y leyenda
ggplot(data, aes(y = Label, x = Value, color = abs(Value - central) < (UpperCI - LowerCI) / 2)) +
  geom_vline(xintercept = valor_central, linetype = "dotted", color = "black", size = 1.5) +
  geom_point(size = 3) +
  scale_color_manual(values = c("FALSE" = "red", "TRUE" = "green"),
                     labels = c("FALSE" = "out", "TRUE" = "in")) +
  geom_segment(aes(xend = UpperCI, yend = Label, x = LowerCI), lineend = "square") +
  
  labs(
    title = "Gráfico de Intervalo de Confianza",
    x = "Valor",
    y = "Etiqueta"
  ) +
  guides(color = guide_legend(title = "Tipo"))  # Título de la leyenda

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

#########5000 Simulaciones#############
# Paso 1: Generar datos de una distribución exponencial (población)
tamano_poblacion <- 10000
poblacion <- rexp(tamano_poblacion, rate = 0.5)
media_poblacion <- mean(poblacion)

# Paso 2: Realizar una sampling distribution de la media
tamano_muestra <- 30
num_simulaciones <- 5000
medias_muestrales <- c()

for (i in 1:num_simulaciones) {
  muestra_indices <- sample(poblacion, size = tamano_muestra, replace = FALSE)
  medias_muestrales[i] <- mean(muestra_indices)
}


# Paso 3: Calcular el intervalo de confianza del 95%
alpha <- 0.05
sigma <- 2
se <- sigma/sqrt(tamano_muestra)
error <- qnorm(1 - alpha/2) * se
left <- medias_muestrales - error
right <- medias_muestrales + error


# Paso 4: Crear un dataframe para ggplot2
data <- data.frame(
  Value = medias_muestrales,
  Label = c(1:num_simulaciones),
  LowerCI = left,
  UpperCI = right,
  central = media_poblacion  # Media de los datos (mu)
)
valor_central = data$central

# Paso 5: Graficar las medias obtenidas junto con sus intervalos de confianza
# Crea el gráfico de intervalo de confianza en vertical con segmentos de error y leyenda
ggplot(data, aes(y = Label, x = Value, color = abs(Value - central) < (UpperCI - LowerCI) / 2)) +
  geom_vline(xintercept = valor_central, linetype = "dotted", color = "black", size = 1.5) +
  geom_point(size = 3) +
  scale_color_manual(values = c("FALSE" = "red", "TRUE" = "green"),
                     labels = c("FALSE" = "out", "TRUE" = "in")) +
  geom_segment(aes(xend = UpperCI, yend = Label, x = LowerCI), lineend = "square") +
  
  labs(
    title = "Gráfico de Intervalo de Confianza",
    x = "Valor",
    y = "Etiqueta"
  ) +
  guides(color = guide_legend(title = "Tipo"))  # Título de la leyenda

# Contar los valores en verde (dentro del intervalo de confianza)
valores_verdes <- sum(abs(data$Value - data$central) < (data$UpperCI - data$LowerCI) / 2)

# Calcular la proporción de valores en verde
proporcion_verdes <- valores_verdes / num_simulaciones

# Imprimir los resultados
cat("Valores en verde (dentro del intervalo de confianza):", valores_verdes, "\n")

## Valores en verde (dentro del intervalo de confianza): 4749

cat("Proporción de valores en verde:", proporcion_verdes, "\n")

## Proporción de valores en verde: 0.9498

Como se ve en los resultados numericos el porcentaje de valores en verde es muy similar al valor del intervalo de confianza solicitado, y cada vez que se repite el experimento se acerca siempre.

Pregunta 3: Estimación de Máxima Verosimilitud

En esta pregunta deberán trabajar con el dataset Body Measurements_original.csv. El objetivo será visualizar e inferir los parámetros que componen a dos variables del dataset. Para esto deberá visualizar el comportamiento de la likelihood, utilizando diferentes cantidades de datos, y realizar la optimización de la likelihood para obtener los estimadores de las variables a través de la función nlminb(). Notar que esta pregunta consiste en dos partes.

1. La primera actividad consiste en realizar inferencia sobre la variable TotalHeight, donde deberá asumir y realizar los siguientes puntos:

• Asuma que los datos de TotalHeight distribuyen de forma gaussiana. Visualice esto a través de un gráfico de densidad en la variable TotalHeight.
• Grafique a través de un gráfico de calor el rango de valores en que se mueve la solución del problema de likelihood, para esto defina su problema en base a -log(likelihood). Mas información se entrega en el esqueleto de la pregunta.
• Visualice cómo se comporta la -log(likelihood) con diferentes cantidades de datos para la estimación de \(\mu\); para visualizar esto, fije \(\sigma\) en 12 y varié solamente el valor de \(\mu\). Para observar el comportamiento del estimador con diferentes cantidades de datos, realice un subsampling de 100 datos, 300 datos y la totalidad de los datos de la variable TotalHeight. ¿Qué se puede observar acerca del comportamiento del estimador \(\mu\)?, responda brevemente esta pregunta.
• Finalmente aplique el comando nlminb() sobre la likelihood y encuentre el máximo o mínimo del problema a optimizar. Señale implícitamente cuales son los valores obtenidos para \(\sigma\) y \(\mu\).

2. Para la segunda parte deberá escoger otra de las variables que componen el dataset, como por ejemplo Age, y estimar a traves de la -log(likelihood) solo los parámetros de la distribución que observa (notar que solo debe inferir los estimadores de variable escogida). Para señalar la distribución de los datos se recomienda realizar el plot de densidad y comparar con el comportamiento de las distribuciones teóricas vistas en clases.

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo, la librería MASS).

Respuesta

Primero, visualicemos la densidad de la variable TotalHeight

data_raw <- read.csv("Body Measurements _ original.csv")

totalheight_plot <- ggplot(
  data = data_raw,
  aes(x = TotalHeight, color = "density")
)

totalheight_dens <- geom_density(color = "#67B7D1")

totalheight_hist <- geom_histogram(
  aes(y = after_stat(density)),
  bins = 10,
  fill = "#67B7D1",
  alpha = 0.5
)

final_plot <- totalheight_plot + totalheight_dens + totalheight_hist +
  ylab("") + xlab("") + scale_color_manual(values = c("density" = "#67B7D1"))


ggplotly(final_plot) %>%
  layout(
    xaxis = list(title = "TotalHeight"),
    yaxis = list(title = "Frequency")
  )

Claramente se aprecia que esta variable se asemeja mucho a una distribución normal. Ahora, veamos un gráfico de calor para visualizar el comportamiento de \(\mu\) y \(\sigma\).

likelihood <- function(mu, sigma, sub = NULL) {
  x_org <- data_raw$TotalHeight
  if (is.null(sub))
    x <- x_org
  else
    x <- sample(x_org, sub)
  n <- length(x)
  x_mean <- mean(x)
  x_var <- var(x)
  first_sum <- (n * x_var) / (2 * sigma^2)
  second_sum <- (n * (x_mean - mu)^2) / (2 * sigma^2)

  -n * log(sigma) - first_sum - second_sum
}

mu <- seq(from = 20, to = 80, by = 0.5)
sigma <- seq(from = 5, to = 23, by = 0.5)
ll_plot <- Vectorize(likelihood)
ll_plot <- outer(X = mu, Y = sigma, ll_plot)

filled.contour(x = mu, y = sigma, z = ll_plot,
               xlab = expression(mu), ylab = expression(sigma))

Se puede apreciar que, muy probablemente, se tiene \(\mu \in (40, 60)\) y \(\sigma \in (5, 10)\).

Ahora, fijemos \(\sigma = 12\) y variemos el valor de \(\mu\)

sigma <- 12
p <- plot_ly(x = mu, y = likelihood(mu, sigma, 100),
             type = "scatter", name = "100 datos")
p <- add_trace(p, y = likelihood(mu, sigma, 300), name = "300 datos")
p <- add_trace(p, y = likelihood(mu, sigma), name = "Todos los datos")
p <- p %>% layout(title = "Likelihood para valores de mu y sigma = 12",
  xaxis = list(title = "mu"),
  yaxis = list(title = "valor de likelihood")
)
p

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

Se puede apreciar que el máximo valor de la función likelihood disminuye mientras aumenta la cantidad de datos considerada. No obstante, esto no influye en el óptimo para \(\mu\), ya que este sigue siendo el mismo para cualquier cantidad de datos ya que la variable es la misma.

Minimicemos la función likelihood para encontrar los valores de \(\mu\) y \(\sigma\).

# Primero, redefinimos la función likelihood para aplicar nlminb
likelihood <- function(param) {
  mu <- param[1]
  sigma <- param[2]
  x <- data_raw$TotalHeight
  n <- length(x)
  x_mean <- mean(x)
  x_var <- var(x)
  first_sum <- (n * x_var) / (2 * sigma^2)
  second_sum <- (n * (x_mean - mu)^2) / (2 * sigma^2)

  n * log(sigma) + first_sum + second_sum
}

# Aplicamos nlminb
MLE <- nlminb(
  objective = likelihood, 
  start = c(10, 10), 
  lower = c(0.001, 0.001), 
  upper = c(500, 500)
)
cat("mu = ", MLE$par[1], ", sigma = ", MLE$par[2])

## mu =  48.1187 , sigma =  12.15672

De aquí obtenemos que los valores son \(\mu = 48.118\) y \(\sigma = 12.156\).

Ahora, trabajaremos con la variable ShoulderToWaist. Veamos su densidad:

data_plot <- ggplot(
  data = data_raw,
  aes(x = ShoulderToWaist, color = "density")
)

data_dens <- geom_density(color = "#67B7D1")

data_hist <- geom_histogram(
  aes(y = after_stat(density)),
  bins = 10,
  fill = "#67B7D1",
  alpha = 0.5
)

final_plot <- data_plot + data_dens + data_hist +
  ylab("") + xlab("") + scale_color_manual(values = c("density" = "#67B7D1"))


ggplotly(final_plot) %>%
  layout(
    xaxis = list(title = "ShoulderToWaist"),
    yaxis = list(title = "Frequency")
  )

Notemos que esta variable también se asemeja a una distribución normal. Usando el mismo método, se ve el mapa de calor y se calculan los valores óptimos.

likelihood <- function(mu, sigma) {
  x <- data_raw$ShoulderToWaist
  sum(log(dnorm(x, mean = mu, sd = sigma)))
}

mu <- seq(from = 10, to = 30, by = 0.5)
sigma <- seq(from = 3, to = 10, by = 0.5)
ll_plot <- Vectorize(likelihood)
ll_plot <- outer(X = mu, Y = sigma, ll_plot)

filled.contour(x = mu, y = sigma, z = ll_plot,
               xlab = expression(mu), ylab = expression(sigma))

# Se redefine la función likelihood para usar nlminb
likelihood <- function(param) {
  mu <- param[1]
  sigma <- param[2]
  x <- data_raw$ShoulderToWaist
  -sum(log(dnorm(x, mean = mu, sd = sigma)))
}

MLE <- nlminb(objective = likelihood, start = c(10, 10), lower = c(0.001, 0.001), upper = c(500, 500))

cat("mu = ", MLE$par[1], ", sigma = ", MLE$par[2])

## mu =  17.90084 , sigma =  5.375555

De aquí obtenemos que los valores encontrados son \(\mu = 17.9\) y \(\sigma = 5.375\).

---

Pregunta 4: Bootstrap

Para esta pregunta será necesario cargar el dataset SAT_GPA.csv, con el que estudiaremos la correlación entre las variables SAT y GPA. Dentro de las variables: GPA representa el rendimiento académico de un estudiante en el sistema estadounidense, mientras que SAT es una prueba de admisión universitaria en estados unidos. Las actividades por realizar en esta pregunta son las siguientes:

• Visualizar la correlación de los datos.
• Obtener la correlación entre los datos de las variables Total, quien representa el resultado de la prueba SAT en USA, y la variable GPA.
• Programar el método Bootstrap para calcular el error estándar de la correlación. Para la simulación utilice un numero de muestras (B) igual a 5000 o algún otro numero de este orden.
• Visualizar a través de un gráfico el histograma obtenido al realizar el muestreo con Bootstrap.
• Obtener el 95% de intervalo de confianza de la estimación por cuantiles como se vio en clases.

Nota: No se permite la utilización de librerías de bootstrap para esta parte.

Respuesta:

my.frame <- read.table(file = "SAT_GPA.csv",header = TRUE,sep = ",", fileEncoding = "latin1")
# Obtener los valores de las variables SAT y GPA
sat_scores <- my.frame$Total
gpa_scores <- my.frame$Gpa

# 1. Visualizar la correlación de los datos
correlation_data <- cor(my.frame)
print(correlation_data)

##                 Rank Bangla     English     Botany    Zoology     Physics
## Rank       1.0000000     NA -0.37957891 -0.6617171 -0.5187892 -0.62425835
## Bangla            NA      1          NA         NA         NA          NA
## English   -0.3795789     NA  1.00000000  0.2034019  0.2517617  0.09586337
## Botany    -0.6617171     NA  0.20340190  1.0000000  0.4981838  0.33292199
## Zoology   -0.5187892     NA  0.25176174  0.4981838  1.0000000  0.22000614
## Physics   -0.6242584     NA  0.09586337  0.3329220  0.2200061  1.00000000
## Chemistry -0.5963821     NA  0.25630775  0.2648451  0.2740580  0.40822483
## Math      -0.5194478     NA  0.06797163  0.1648513  0.1381620  0.31918558
## Ict       -0.5183436     NA  0.28033712  0.3222062  0.1219240  0.22727532
## Total     -0.9444559     NA  0.46241917  0.6996960  0.5443674  0.62214193
## Gpa       -0.9810238     NA  0.39771083  0.6684849  0.5228122  0.62404239
##            Chemistry        Math        Ict      Total        Gpa
## Rank      -0.5963821 -0.51944782 -0.5183436 -0.9444559 -0.9810238
## Bangla            NA          NA         NA         NA         NA
## English    0.2563078  0.06797163  0.2803371  0.4624192  0.3977108
## Botany     0.2648451  0.16485135  0.3222062  0.6996960  0.6684849
## Zoology    0.2740580  0.13816205  0.1219240  0.5443674  0.5228122
## Physics    0.4082248  0.31918558  0.2272753  0.6221419  0.6240424
## Chemistry  1.0000000  0.26818598  0.2336958  0.5778151  0.5910724
## Math       0.2681860  1.00000000  0.1674076  0.4989442  0.5559980
## Ict        0.2336958  0.16740757  1.0000000  0.5455583  0.5146484
## Total      0.5778151  0.49894420  0.5455583  1.0000000  0.9654575
## Gpa        0.5910724  0.55599797  0.5146484  0.9654575  1.0000000

cat("\n")

# 2. Visualizar la correlación de los datos
correlation_var <- cor(sat_scores, gpa_scores)
cat("Correlación entre SAT y GPA:", correlation_var, "\n")

## Correlación entre SAT y GPA: 0.9654575

# 3. Metodo bootstrap para calcular el error estandar de la correlacion
myboot<-function(x,y,fun,nRuns,sampleSize,alpha){
  
  values<-vector()
  
  for(i in 1:nRuns){
    
    samp.i<-sample(x, size = sampleSize, replace = TRUE)
    samp.j<-sample(y, size = sampleSize, replace = TRUE)
    values[i]<-fun(samp.i,samp.j)
    
  }
  

  point.est <-fun(x,y)
  
  # error estandar
  se <- sd(values)
  
  # 4. Visualizar el grafico de Histograma
  hist(values,main = "Histograma de correlaciones", xlab="Correlacion",ylab = "Frecuencia")
  
  
  #5. Obtener el 95% de intervalo de confianza de la estimación por cuantiles
  l.CI <- quantile(values, alpha/2)
  
  u.CI <- quantile(values, 1-alpha/2)

  
  return(c("Point Estimate"=point.est,
  "Standard error"=se,
  "Lower CI limit" = l.CI,
  "Upper CI limit" = u.CI))
}

cat("\n")

myboot(sat_scores,gpa_scores,cor,5000,100,0.05)

##       Point Estimate       Standard error  Lower CI limit.2.5% 
##            0.9654575            0.1002537           -0.1952264 
## Upper CI limit.97.5% 
##            0.1966485

---

 

---

A work by CC6104